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Abstract 

(N 

A geometrically exact membrane formulation is presented that is based on curvilinear coordi- 
nates and isogeometric finite elements, and is suitable for both solid and liquid membranes. 
The curvilinear coordinate system is used to describe both the theory and the finite element 
equations of the membrane. In the latter case this avoids the use of local cartesian coordinates 
at the element level. Consequently, no transformation of derivatives is required. The formula- 
te tion considers a split of the in-plane and out-of-plane membrane contributions, which allows the 
construction of a stable formulation for liquid membranes with constant surface tension. The 
proposed membrane formulation is general, and accounts for dead and live loading, as well as 
enclosed volume, area, and contact constraints. The new formulation is illustrated by several 
challenging examples, considering linear and quadratic Lagrange elements, as well as isogeomet- 
Q ric elements based on quadratic NURBS and cubic T-splines. It is seen that the isogeometric 
elements are much more accurate than standard Lagrange elements. The gain is especially large 
y—^ for the liquid membrane formulation since it depends explicitly on the surface curvature. 
> 

T— I Keywords: contact constraints, curvilinear coordinates, isogeometric analysis, nonlinear finite 

0^ element methods, follower loads, volume constraints. 

d 

^ 1 Introduction 

.!_, Membranes are computationally challenging structures. Their geometry can be complex, they 

^ may undergo large deformations, large rotations and large strains - and thereby behave highly 

nonlinear - and they are characterized by several physical instabilities: They are unstable in 
compression, unstable for out-of-plane loading (in the case of zero in-plane tension) , unstable for 
pressure loading (in the case of rubber membranes) and unstable w.r.t. in-plane loading (in the 
case of liquid membranes). The aim of this paper is to formulate a general, 3D, geometrically 
exact and fully nonlinear membrane model that accounts for pressure loading as well as volume, 
area, and contact constraints and is suitable for both solid and liquid membranes. Our focus is 
on pure membranes, i.e. curved, surface structures that do not support in-plane compression, 
out-of-plane bending, and shearj^ Such a restricted focus is useful due to the large range of 
membrane applications: they appear as inflatable and pressurized structures, like balloons, tubes 
and airbags; as fabrics, tents, canopies, parachutes, foils and sails; as water-filled membrane 

^corresponding author, email: sauer@aices.rwth-aaclien.de 

^We note that in the literature, the term membrane is often also used for the special case of 2D plane-stress 
structures. 
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structures, like inflatable dams; as biological membranes, like blood vessels, cell, diaphragms, 
aneurysms and lung alveoli; as liquid droplets, menisci, bubbles, foams and sprays; as thin sheets 
and films - both liquid and solid - as atomistic membranes, like graphene sheets; as interacting 
membranes, e.g. adhering cells; and in the topic of form-finding and minimal surfaces. 

Computational formulations for 3D, nonlinear membrane go back to the seminal work of Oden 



(Oden and Sate 


) (1967 


), see also 


Oden (: 


lOOi 




Since then, the field has been continuously 


advanced, among others by Fried 


(1982) 


Tang 


(1982 


); 


Roddeman et al. 


(1987); 


Contri and 


Schrefler (1988) 


Wriggers and Taylor ( 199 


0);] 


brahimbegovic and Gruttmann 


(1993 




iaseganu 


and Steigmann 


1994); 


Gosling and Lewis 


(1996) 


Muttin 


(1996); Wu et al. 


(1 


996');'Bonet et al. 


f2000|) ; ,Rumpel and Schweizerhof 


(2003); 


Stanuszek ( 


2003); Weinberg and Neff, (200^ 


i). Many 



of these works are concerned with the topic of wrinkling due to in-plane compression. More 
recently, computational formulations based on curvilinear coordinates have been considered rig- 



orously, both for membranes ( Ambroziak and Klosowski 2006 ) and shells ( Arciniega and Reddy 



2007). Another recent development are rotation-free shell formulations, as they have been con- 



sidered by Flores and Estrada (2007); Linhard et al. (2007); Dung and Wells (2008) and recently 



Benson et al. (2011); Nguyen-Thanh et al. (2011) for isogeometric analysis. Isogeometric for- 
mulations allow the formulation of C^-continuous surface formulations that are advantageous 



for flow simulations (Kiendl et al. , 2010) and sliding contact (De Lorenzis et al. , 2011 Temizer 



et al. , 2012), see also Sauer (2011, 2012) for Hermite-based, C^-continuous contact surfaces 



Relevant to membranes is also the topic of live pressure loading (Buffer, 1984 Schweizerhof and 



Ramm, 1984). Membranes are also an interesting subject in shape optimization (Bletzinger 



et al. 2005 Manh et al. 2011) 



The presented formulation contains several merits and novelties: It allows a split between in- 
plane and out-of-plane contributions, which is used to construct a new formulation for liquid 
membranes. It admits arbitrary elastic material models for solid and liquid membranes. It 
is based purely on displacement-based finite elements and can be used with any kind of such 
elements. It includes, in particular, isogeometric NURBS elements to capture the deforming 
surface geometry to high- accuracy, even for comparably coarse discretizations. It is straight 
forward to implement in an existing FE framework. It avoids the need to transform derivatives 
between configurations and avoids the use of local cartesian coordinate systems. Shells models 
are often formulated using local cartesian coordinate systems, as this allows using classical 



constitutive relations formulated in this manner (Wriggers, 2008). To our mind, there is no 
need for such a detour: The balance laws, kinematics, constitutive relations as well as the FE 
weak forms and corresponding FE arrays can all be formulated efficiently in the curvilinear 
coordinate system. The capabilities of the presented formulation are demonstrated by several 
challenging computational examples, considering pressure loading, inflation and contact. 

The following section presents the theory of nonlinear membranes in the framework of curvilinear 
coordinates, considering pressure loading, volume, area and contact constraints. Sec. [3] proposes 
a straight-forward finite element implementation of the theory, and Sec. |4] presents several 
examples of solid and liquid membranes to illustrate the capabilities of the present formulation. 



2 Nonlinear membranes 

In this section, we summarize the theory of nonlinear membranes in the framework of curvilinear 
coordinates. The membrane kinematics, constitution and balance laws in strong and weak form 
are discussed, and various kinds of constraints are considered. 
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2.1 Surface description in curvilinear coordinates 



The membrane surface, denoted S, is fully characterized by the parametric description 

x = x{e,e). (1) 

This corresponds to a mapping of the point (C^, C^) in the parameter domain V to the material 
point X £ S. In the following, Greek letters are used to denote the two indices 1 and 2. 
Summation is then implied on repeated indices. The tangent vectors to coordinate at point 
X £ S are given by 

dx 

a. = ^, a = l,2 (2) 

The two vectors form a basis for the tangent plane of S. In general, they are not orthonormal. 
This apparent drawback of the description is actually an advantage when it comes to the 
kinematical description. This turns out to be very straightforward, e.g. see Eq. (16). The 
basis at x is characterized by the metric tensor, that has the co-variant components 

aai3 ■= . (3) 

From the inversion 

K^] := (4) 

we obtain the contra-variant components of the metric tensor. Explicitly, we have a^^ = 
aii/detoQ^, a}'^ = — ai2/detaQ^ and a?"^ = 022 / det a^/j . With these a dual basis can be 
constructed: From the co-variant base vectors the contra-variant counterparts 

a" := a"^a^ (5) 

can be determined. We note that summation is implied on repeated indices. Note that a° -a^ = 
a"^ and a^- = ^a, where 5a is the Kronecker symbol. The unit normal of 5 at a; is given by 

ai X a2 

^ = n 11 • (6) 

||ai X a2\\ 

It can be shown that 

||ai X a2\\ = det aap ■ (7) 
The bases {ai, 02, n} and {a^, a?, n} can then be used to decompose any vector v on 5, i.e. 

V = v"' aa + VriU = Vaa"" + v-a_n , (8) 

where Va denote the co- variant, and f " the contra- variant components of v. The derivative of 
the tangent vectors is given by 

daa , , 

Further, we require the so-called co- variant derivative of a^, which is defined by 

:= - r^^a7 (10) 

where F^^ are the Christoffel symbols of the second kind given by F^^ = a^^p ■ a? ■ Introducing 
the identity tensor on S 

\ = aa®a°' = a°'®aa = ^ — n^n, (11) 
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where i is the usual identity tensor in , |^ we can write 

aa;i3 = (n (g) n) a^,^ . (12) 

Contracting with n then yields 

n ■ a^-p = n ■ a^^p = h^p (13) 

which are the co- variant components of the curvature tensor h = a"^ a^. The eigenvalues 
of this tensor are the principal curvatures of surface S. 



2.2 Membrane kinematics 

Next, we consider the deformation of the membrane surface. We therefore distinguish between 
the deformed, current configuration S and the undeformed, initial configuration Sq, see Fig. [T] 



Both surfaces are described by the relations of Sec. 2.1 For surface S we use the lower case 




Figure 1: Mapping between parameter domain V, reference surface Sq and current surface S 



symbols x, a^, aap, a", n and 6^/3 • For surface So we use the corresponding upper case symbols 
Aq,, Aq,/^, and -^|^ In order to characterize the deformation between surfaces 5o and S 
consider the line element 

da;=— dr = a„dr (14) 
and likewise dX = d^" . Contracting with yields d^" = A" • dX , so that 

da; = (a„ ® A") dX . (15) 

Here the tensor 

F = a„ A" (16) 

tilde is used here to indicate standard tensors in 
''Here, the curvature tensor [fca^g] is only needed on S. 
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is the surface deformation gradient of the mapping X ^ x. Likewise we find F ^ 
Through F we thus have the fohowing transformations 



a" 



FAa , 



-T Aa 



Act — F CLd 1 



F'^a" 



Aa a". 



(17) 



Given F, we can introduce the right and left Cauchy-Green surface tensors and their inverses, 
C = F^F = ar.RA''^Al^ , C'^ = a'^^ A^ Af. 



B = FF^ 



■ 



a 



B 



-1 



a 



(18) 



Next we discuss the surface stretch between surfaces Sq and S. The area element da C 5 is 
defined by 



da 



lai X a2\\ dD 



(19) 



{aide) X {a2de)\\ 

where dD := d^^d^^. A corresponding statement follows for dA C 5o- In view of Eq. ([t]) we 
thus have the relations 

det Af^p , 
y/ det aai3 , 



d^l 


= JAdn , 


Ja 


da 


= Jadn, 


Ja 


da 


= JdA , 


J 



(20) 



Ja/JA ■ 



2.3 Momentum balance for membranes 



From the balance of linear momentum the strong form equilibrium equation 

tZ + f = 0, 



(21) 



at a; G 5 can be obtained (Steigmann 1999). Here / is a distributed surface force, that can be 
decomposed as 

f = faa'' +pn = ra^+pn (22) 

where fa and /" are the co-variant and contra-variant in-plane components of / and p is 
the out-of-plane pressure acting on S. Further, denotes the internal traction acting on the 
internal surface _L a". According to Cauchy's formula 



cra^ 



(23) 



where a denotes the Cauchy stress tensor at x £ S, which we consider to be symmetric. We 
emphasize that is not a physical traction since a" is usually not normalized. In general, the 
stress tensor takes the form 



ap + a {n ^ aa + da ^ n) + a n®n 



(24) 



For membranes it is typically assumed that cr'^" = a^^ = 0, so that 



cr = cr"^ aa 



ar- 



In this case we find that = a^"ap such that the co- variant derivative of becomes 

t'^ = a%ap + aP^bapn (25) 



according to Eq. (13). Equilibrium equation (21) thus decomposes into 



fj^fa + = , (in-plane equilibrium) , 

bap + P = , (out-of-plane equilibrium) . 



(26) 
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To close the problem, the usual Dh'ichlet and Neumann boundary conditions 

u = u on d.iS 
t = t on c/fO 

are considered on the membrane boundary dS = duS U dtS. Here, we suppose that the 
prescribed traction i = t" is tangent to S, since out-of-plane boundary forces, as well as 
out-of-plane line and point loads within the surface, lead to singularities in the membrane 
deformation and are therefore not considered in the present formulation. The traction on 
boundary dtS, according to Cauchy's formula, is given by 

t = am , (28) 

where m = maO," is the outward unit normal of dtS. It follows that t = mat"'- 



2.4 Membrane constitution 

The known 3D constitutive models can be adapted to the membrane. We therefore suppose 
a general elastic material relation of the form or = cr(S). For membranes it is useful to 
consider the decomposition B = S + Ag (n (gi n), where A3 is the out-of-plane stretch, and 
cr = cr/t + (T33 (n (g) n), where cr is defined as the in-plane stress tensor (with units force per 
length) and t = A3T denotes the current membrane thickness, for a given reference thickness 
T. Out of these considerations, a relation between the membrane quantities B and cr can be 
obtained. As an example we consider an incompressible Neo-Hooke material, given by 

& = fLB + ql, (29) 

where fi is the shear modulus and q denotes the Lagrange multiplier associated with the incom- 
pressibility constraint. For membranes, the model decomposes into 

a = ifiB + ql)T/\s , 
o"33 = /^A§ + q . 

For incompressibility det B = (JAs)^ = 1. Under the plane stress assumption (T33 = 0, we then 
find q = — /i/J^ and consequently 

a = fi/j{B - 1/J2) , (31) 
with fj, := fiT. Componentwise, in the aa basis, this becomes 

= ^/J[A^P - a^P/j'^) . (32) 
Contracting with a^^, the components o"^ and Gap can be obtainedj^ 

Another example are liquid, e.g. water, membranes governed by constant isotropic surface ten- 
sion 7. In that case a constant stress tensor of the form 

CT^ = 7 5^ (33) 



is obtained. It can be seen that the in-plane equilibrium equation (26 1) is only satisfied for 
/" = 0. This implies that static water membranes cannot equilibrate in-plane loads, and 
are therefore unstable in-plane; a property that needs to be addressed in a computational 



formulation (see Sec.ra. The out-of-plane equation (26 2) now yields 



2H-f + p = 0, 2H:=b%, (34) 

which is the well known Young-Laplace equation. A prominent feature of liquid membranes is 
that they form distinct contact angles. This property is not addressed here. 



'Due to the symmetry of a"^ the ordering of indices does not matter in ctS, i.e. <j' 
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2.5 Membrane weak form 



Next, we derive the weak form corresponding to equilibrium equation (21). Consider a kine- 
matically admissible variation of S, denoted w G W, where W denotes a suitable space for w. 

(35) 



Contracting Eq. (21) with w and integrating over S yields 

w ■ {f^^ + /) da = yweW 



L 

Considering w = Wa a" + wn, this expands into 

j + r) da + j w (a"'' 6c./3 +p) da = Vii; G W , 



(36) 



i.e. it splits into the in-plane and out-of-plane parts identified in Eq. (26). Such a split is useful if 
different approximation techniques are chosen for the in- plane and out-of-plane response. Using 
the divergence theorem for curved surfaces (Gurtin and Murdoch, 1975), the first in-plane term 
is rewritten into 



Wa;l3 'y'^^ da 



/ Wq, cj"'^^ da = / (wq, o""^).^ da — / 
J S J S J s 

= I Wa cr°^ m/3 ds — Wa-(5 (7°"^ da , 
JdS Js 



(37) 



where rua = m ■ aa are the co-variant components of the unit normal m on the line dS. Since 
u;" = on duS and since WaC^'^^mi^ = Wa on dtS expression (36) thus becomes 



Gint ~ Gext 



G W 



where 



G 



int 



a-p cj"'^ da — w a'^^ b^p da 



ext 



Wa f"' da + I Wat°'ds+ I wpda , 

dtS 



(38) 



(39) 



are the internal and external virtual work contribution due to variation w. Considering Wa = 
and w = subsequently, the weak form can be decomposed into the weak forms 



Wa-R cr"''^ da - / Wa P da 



Wat"' ds = y Wa ^ y^a (iu-plaue) 



dtS 



/ w (7°^ bai3 da+ wpda = \/wG Wn (out-of plane) 
Js Js 



(40) 



Such a split is advantageous for the description of liquid membranes. Since liquid membranes 
are inherently unstable in-plane, they can be stabilized by providing additional stiffness via 
Eq. (40 1) without affecting the out-of-plane response. This is demonstrated in the examples of 



(41) 



Sec. 133] and I45l 
Otherwise, considering 

Wa-I3 = {W ■ aa)-l3 = W-ji ■ tta + W hap , 

the two terms of Gint can be combined into 



G 



int 



w-ry ■ a'^^ as da . 



(42) 
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It is noted that for this expression only single derivatives of variation w and configuration x are 



required, while in the decomposed formulation of Eq. (40) second derivatives of x appear. Intro- 
ducing the surface Kirchhoff stress tensor r = Jcr, which eliminates one J from expression (32), 
the last equation can be rewritten into 



Gint= / w.^-T^'PapdA . (43) 



Within framework (38), both dead and live loading can be considered. This is discussed further 



in Sec. 3.2 Beforehand, we discuss several useful membrane constraints. 



2.6 Volume constraints 

The volume of the domain D enclosed by the membrane may be constrained. An example is a 
cell containing incompressible fluid. Formally, the volume constraint is written as 

g,:=V-Vo = 0, (44) 

where Vq and V denote the initial and current volumes enclosed by the initial and current 
membrane configurations. Since dv = x-n da/3 and d^ = X ■ N dA/3, these can be computed 
by the surface integration 

V = l[x-nda, Vo = l [ X-NdA (45) 

"3 J5 -3 J So 

These expressions are valid for closed surfaces, and care has to be taken when modeling open 
membranes. In this case one must account for the volume contribution associated with the miss- 



ing surface. Eq. (44 ) can be included in the formulation by the Lagrange multiplier method. The 



Lagrange multiplier associated with the volume constraint is the internal membrane pressure p. 



Remark: We note that the governing equations (38) and (44) can be derived from a variational 



principle for conservative systems. This is the case for the constitutive models discussed in 



Sec. 2.4 and for pressure loading of closed membranes surfaces. 



2.7 Area constraints 

Another useful constraint is a constraint on the membrane surface area. For example, red blood 



cells are known to conserve the surface area during deformation (Kloeppel and Wall, 2011). 
Formally this is expressed as 

g^:=A-Ao = 0, (46) 

with 



A= da, Ao= dA. (47) 

Js JSo 

The area constraint is not considered further in this paper. In principle, they can be treated in 
an equivalent manner to volume constraints. 



2.8 Contact constraints 

Contact is characterized by the impenetrability constraint 

5n < , (48) 
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where 



9n ■■= {x 



(49) 



denotes the normal gap between the membrane point x & S and the surface F of a neighboring 
obstacle. Here, the unit vector rip denotes the surface normal of T at the point Xp, which is the 
solution of the minimum distance problem 



|y| min [x — y) for x G 5| 



(50) 



We note that this minimization can cause difficulties for complex surface geometries (Wriggers 



2006). Constraint (48) can be included in the membrane formulation by various methods. The 



simplest of these is the penalty formulation. In this case the contact traction tc acting at x £ S 
is given by 

-enS'n^-p , ffn < , 
, 5n > , 



(51) 



where en is the penalty parameter. The contact forces contribute to virtual work balance (38). 
This can be expressed by including 



Gr 



w ■ tcda 



(52) 



on the right hand side of Eq. (38). For two deformable membranes in contact, weak form (38) 



must be satisfied for each membrane, and contribution Gc is added correspondingly to each 
weak form. To avoid a surface bias it is advantageous to treat both contact pairs equivalently 
as is done in the two-half-pass algorithm of Sauer and De Lorenzis (2012) . For some problems, 



like adhesion, the contact constraint is replaced by suitable constitutive contact laws of the form 
tc = icign), see Sauer and Li (2007). 



3 Finite element discretization 



The governing equations (38) and (44) are solved by the finite element (FE) method. The ini- 
tial surface Sq is therefore discretized into a set of finite elements that are defined by nodal 
points Xj or control points in the case of isogeometric FE. The deforming membrane is then 
described by the motion of the nodal points Xj — t- xj, which corresponds to a Lagrangian FE 
description. The deformed configuration of element is denoted fi^. Here we consider quadri- 
lateral elements since these can be conveniently related to a master element in the parameter 
domain G [—1, 1]. 



3.1 Finite element interpolation 

Within elements Og geometry is approximated by the nodal interpolations 

X^X^ = Y,N,Xj, (53) 

and 

X ^ x'^ = ^NiXi , (54) 
I 

where Nj = A'^/(^^,^^) denotes the nodal shape function defined on the master element in 
parameter space. The summation is carried out over the n^e nodes of the element. Here, the 
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following quadrilateral elements are considered: 4-noded linear Lagrange elements, 9-noded 
quadratic Lagrange elements, quadratic NURBS elements, and T-spline elements. In principle, 
any other element type can also be considered. For isogeometric elements the shape functions 



are constructed via the Bezier extraction operation (Borden et al. , 2011 Scott et al. 2011) 



According to Eq. ([2]), the tangent vectors are thus approximated by 



XI 



(55) 



where Nj^^ = dNj/d^,". Considering a Bubnov-Galerkin formulation the variation w is approx- 
imated in the same way as the deformation, i.e. 



w 



I 



For shorthand notation, we rewrite Eqs. (53), (54) and (56) 

X ^ NX. 



e , X ^ NXe , W ^ NWe , 



(56) 



(57) 



where N := [A^il, Njl, ...] is a (3 x 3n^e) array with the usual identity tensor 1 and Xe, Xe, 
and We are vectors containing the stacked nodal values for the element]^ In order to discretize 
the weak form we need to discretize Wa, w and w-a- We find 



w 

Wr 



W ■ tta 

w ■ n 



w^N^N^Xe 



(58) 



where N^q := ■■■Nj^al, ...]. Note that for a vector like w, the co-variant derivative 

w-a coincides with the regular partial derivative w^a- The surface normal n is given through 
definition ([6]) and approximation (55). According to (13), the components of the curvature 
tensor become 

ba/S~n- N^al3 Xe . (59) 



With the above expressions, we further find 



Wa;l3 



(^NJ N,„ + N'^ (n n) N^^ib) 



(60) 



according to eq. (41). 



3.2 Discretized weak form 



The above expressions are now used to discretize the membrane weak form of Sec. 2.5 The 
surface integration is carried out over the element domains fi*^ and then summed over all riei 
FE as 



G 



mt — / , '-^int ) 
e=l 



Eg 

e=l 



ext 



For the internal virtual work of eq. (|43|) we now have 



int 



W; 



(61) 



(62) 



^Non-italic discrete arrays X, x, w and N should not be confused with italic field variables X, x, w and AT. 
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according to approximations (55) and (58). Writing Gf^^^ = wj"^^^., we identify the internal FE 
force vector 



(63) 



As noted in Eq. (39 1), the internal virtual work can be split into in-plane and out-of plane 



contributions. At the element level these are 



"-^inti 



e 

into 



(64) 



such that G 



int 



^fnti + ^fnto- -'^^ view of Eqs. (58), (59) and (60), the corresponding force 



vectors become 



^nti = / r"'' (n^ N,;3 + (n ^ n) N,„^) d^l Xe , 
ffnto = - / ^"'^ N^(n ® n) N,„;3 d^ x^ . 



(65) 



For Gext we consider external loading of the form / = /q/J + pn, where /g and p are given 
loading parameters. The first corresponds to a dead force per reference area, the second to a 
live pressure. According to Eq. (39 2) we then have 



"-'ext 



w ■ /q d^ +/ w ■ tds + wpda , 



which yields the external force vector 



N^tds+ / N'^pnda . 



The original weak form (|38|) now yields the descretized version 

, 



w 



'[fi 



int '■cxt 



(66) 



(67) 



(68) 



where fint and fext are obtained from the assembly of the corresponding elemental force vectors, 
and w is the kinematically admissible set of all nodal variations. These are zero for the nodes 
on the Dirichlet boundary duS. For the remaining nodes, Eq. (68) implies 

f • — fint fext 







(69) 



which is the discretized equilibrium equation that needs to be solved for the unknown nodal 
positions x; see Sec. |3.5[ 



We note, that in this formulation no mapping of derivatives between master and current config- 
uration is required. Also no introduction of local, cartesian bases are needed. The formulation 
thus is straight forward and efficient to implement. 



3.3 Contact contributions 



The proposed membrane model can be easily extended to include contact, provided a 3D contact 



algorithm is available. The contact contribution (52) simply yields the force vector 



tc da 



(70) 



that needs to be included in Eq. mQv. For details on the the FE implementation of Eq. (70) we 



refer to Sauer and De Lorenzis ([2012|). 
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3.4 Discretized volume constraint 



The volume, enclosed by the discretized membrane surface, is obtained as 

3 



V 



-V / n^NdaXe 



(71) 



according to Eq. (45). For the volume constraint, = V — Vq = Vq can be considered as an 



externally prescribed volume, e.g. during inflation, or as the initial value of V . 



3.5 Solution method 



The volume constraint is included in the formulation by the Lagrange multiplier method. The 
Lagrange multiplier associated with the constraint is the pressure, p, acting on membrane. 



Combining (44) with (69) leads to the system 



f(x,p) 



, 

, 



(72) 



that needs to be solved for the unknown nodal position x and pressure p. Due to the nonlinear- 
ities of the model, this is solved with Newton's method. Therefore, the linearization of f and 
5v w.r.t. X and p are needed. This is discussed in Appendix [B) 



3.6 Hydrostatic pressure 

In some applications, the pressure p may vary locally. In static examples this is typically due 
to gravity. An example is the hydrostatic pressure distribution in a water-filled membrane. In 
this case, we have 

p = Pv+Ph , (73) 

where is the hydrostatic, height dependent, pressure and Pv is the pressure associated with 
the volume constraint. The former is simply written as 

Ph = -pg-x (74) 

where p is the density of the pressure causing medium, and g is the gravity vector]^ The value 
of Pv is then the (constant) datum pressure at the origin. 



3.7 Numerical quadrature 



In parameter space, each element is defined on the master domain ^" G [—1, 1], a = 1,2. The 
integrals from above are mapped to the master domain using transformations ( 20 ) . Integration 
is then carried out with standard Gaussian quadrature on the master domain. 

^typically g — —[0, 0, , where g is the gravity constant 
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3.8 Monitoring compression 



Membranes do not support in-plane compression. The absence of physical bending stiffness leads 
to buckling of the structure, known as wrinkling in the case of membranes. To avoid membrane 
compression in our formulation during computations, we simply monitor the minimum principal 
stress 

<7min ~ Y ~ V T ~ 

where Ii = tr <t = cr^ and I2 = det cr = det are the two invariants of the surface stress tensor. 
We note that fTmin does not imply the automatic failure of the discretize membrane structure 
as some numerical bending stiffness may be present. More involved wrinkling criteria can be 
found in the literature, seelLu et al.l (|200ip and lYoun and Leel (|2006p . 



(75) 



4 Numerical examples 



The proposed membrane model is illustrated by several examples, considering both solid and 
liquid membranes under dead, pressure, and volume loading. Standard linear and quadratic La- 
grange finite elements as well as quadratic NURBS and cubic T-spline finite elements, providing 
C^- and C^-continuous surface descriptions, respectively, are used. 



4.1 Inflation of a spherical balloon 



We first consider the inflation of a spherical rubber balloon and use it for validation, since an 
analytical solution exists for this problem. The rubber behavior is described by the incompress- 
ible Neo-Hookean material model (29). The finite element model of the balloon, modeled as 
1 /8th of a sphere, is shown in Fig. [2^. Appropriate boundary conditions are provided to main- 




Figure 2: Inflated balloon: (a) initial and current configuration (for V = 10 Vq); (b) pressure- 
volume relation for y G [1 10] Vq (FE result for 3 quadratic FE). 



tain the symmetry of the inflating structure. The relation between current and initial radius is 
denoted r = XR. The circumference of the balloon, proportional to r, is thus stretched by A 
such that the surface deformation gradient is = Al and the area change is given by J = A^. 
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Due to incompressibility this results in the thickness change t = T/J. According to Eq. (31), 
the in-plane normal stress within the balloon thus is a = a/t = /iT(l — X^)/t, which is equal to 
the well-known formula a = pr/2/t. We thus obtain the pressure-stretch relation 

^fi-iV ,76) 



or, equivalently, the pressure-volume relation 

where Vq = 4tt/3R^ is the initial balloon volume. The p — V relation is shown in Fig. ^p. 
The pressure increases quickly, peaks and then decreases gradually. This behavior is typical 
for the inflation of rubber membranes. The FE computation of such problems should therefore 
be carried out by prescribed volume loading instead of prescribed pressure loading. The pro- 
posed FE formulation can capture the analytical behavior very nicely. This is shown by the 
convergence plot of Fig. Isk. Here, the number of Gaussian quadrature points per elements are 




Figure 3: Inflated balloon: convergence of the pressure at V = 10 Vq: (a) convergence with 
mesh size; (b) NURBS convergence with quadrature accuracy. 



2x2 for linear Lagrange and 3x3 for quadratic Lagrange and NURBS elements. Since NURBS 
elements describe the spherical geometry exactly, they can solve the problem exactly with only 
one element - provided sufficiently many quadrature points are used. This is shown in Fig. ^p. 
The results shown here validate the proposed membrane formulation. 

Fig |4] shows the deformed FE meshes and the error in the membrane stress a for the three 
different element types considered here. As is seen, the error is smallest for NURBS FE. 



4.2 Inflation of a square sheet 

As a second example we consider a square membrane sheet with dimension 4Lo x 4Lo, apply 
an isotropic pre-stretch of Aq = 1.05 to provide initial out-of-plane stiffness, and then inflate 
the structure by a prescribed volume, as is shown in Fig. [Sj 

Fig. [6^ shows the pressure-volume relation for the three considered elements. The accuracy is 
highest for NURBS elements and lowest for linear elements. This is seen by the convergence 
behavior of the different element types, shown in Fig. [HJd. 
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Figure 4: Inflated balloon: Error in the in-plane stress a = pr/2 for: (a) 12x8 linear FE, (b) 
3x8 quadratic FE, (c) 1x8 NURBS FE. 
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Figure 5: Inflated square sheet: configurations for V = {0, 1, 2, 3, 4, 6, 8, 10}Vo, where Vq 
4Lq. The coloring shows the area stretch J (which is identical to the thickness decrease). 



Fig. [7] shows the deformed sheet for a prescribed volume oi V = 5000 Vq for the three element 
types. As seen, all element types can accommodate enormous deformations, even for relatively 
coarse meshes. The comparison with the fine NURBS mesh in Fig. [8^ shows that there are 
still considerable inaccuracies present in all three formulations. The NURBS result is fully 
C^-continuous. In the example, particularly large deformation occur at the bottom and in the 
corners of the sheet, as is seen in the close-up of Fig. [Sja. The deformation in the corner shown a 
tendency towards wrinkling. We observed that a further mesh refinement led to non- convergent 
Newton behavior, indicating instabilities. A computational scheme for wrinkling is required to 
handle this case. 



4.3 Contact between balloon and cushion 



The third example considers a spherical, water-filled balloon in contact with a cushion. The 
balloon is loaded by hydrostatic pressure loading. The cushion is modeled by a square sheet that 
is fixed along the boundary and supported by internal pressure arising from constraining the 
volume beneath the sheet. The initial size of the sheet is 2R x 2R, where R is the undeformed 



radius of the balloon. Both, balloon and sheet are modeled by material law (29) considering 
equal fi. They are both pre-stretched isotropically by Aq = 1-1, i.e. the constrained balloon 
volume is Vq = 37r (Ai?)^/4. The problem is computed by gradually increasing the gravity level, 
g, pulling on the water inside the balloon. Quadratic finite elements are used. Contact is 



modeled by the two-half-pass contact algorithm (Sauer and De Lorenzis, 2012) considering the 
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Figure 6: Inflated square sheet: (a) pressure-volume relation; (b) pressure convergence at V = 
5000 Vb (compared to a quadratic NURBS mesh with 23233 dofs). 




Figure 7: Inflated square sheet: deformation at ^ = 5000 Vq for (a) 8 x 8 linear elements, (b) 
4x4 quadratic elements, and (c) 8 x 8 NURBS elements. The color shows the area stretch J 
displayed as log^o J- 



augmented Lagrange multiplier method. In principle, any 3D contact algorithm can be applied 
straight forwardly to the proposed membrane formulation. Fig. [9] shows the deformation of 
balloon and cushion for various gravity levels. As shown, the deformation becomes very large, 
which makes the problem very challenging. The example is interesting as it involves large 
deformations, contact, pressure loading, hydrostatic loading and two volume constraints. 



4.4 Growth of a hemispherical water droplet 



As a validation of the formulation for liquid membranes, we consider the growth of a hemi- 
spherical droplet resting on a rigid substrate and maintaining a contact angle of 90°. The 
problem is similar to the balloon inflation example (Sec. 4.1), and the same FE meshes are 
used. For a liquid water membrane the membrane stress is given by Eq. (33), i.e. the stress is 
deformation independent. This implies that only the out-of-plane but not the in-plane forces 
provide stiffness and the formulation is unstable in itself. The formulation can be stabilized 
by adding deformation dependent in-plane forces through Eq. (65 1). These forces should not 
influence the out-of-plane behavior such that the original liquid membrane formulation remains 
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Figure 8: Inflated square sheet: deformation at V = 5000 Vb for 88 x 88 NURBS elements: (a) 
overall deformation, (b) deformation at the corner. 



unaffected. We simply use the incompressible Neo-Hookean model to provide the additional 
in-plane stability. The Neo-Hookean material parameter then becomes a numerical stability 
parameter that should not affect the physical results. The internal forces acting on the finite 
element nodes are then simply given by 



'■mt — ^mt I '^liquid 



al3 \ 
inti l^solidy 



(78) 



For this example, like in Sec. 4.1 , the pressure- volume relation is also know analytically. Setting 
a = pr/2 = 7, with r = XR and V = X^Vq, we find 



pR 

7 



(79) 



Since the pressure remains positive (and is thus stabilizing the structure) the volume can also 
be decreased. The computed pressure and the convergence of the proposed finite element for- 
mulation to the analytical result are shown in Fig. 10 Several values for the numerical stability 
parameter fj, are considered. They all converge to the desired analytical result. Quadratic 
Lagrange and NURBS elements are considered, and it is seen that the NURBS formulation 
converges much faster. This is attributed to the higher surface continuity that appears in fj^^j 



according to Eq. (65 ,1). Decreasing /i improves the accuracy. A more detailed analysis of the 



model proposed in Eq ( 78 ) along with the effect of parameter /i is required and will be considered 
in the future. 



4.5 Liquid droplet on a rigid substrate 

The last example examines a static water droplet in contact with a rigid substrate. A distinct 
feature of liquid membranes is that they can form sharp contact angles at the contact bound- 
ary. In this case, the membrane surface forms a kink at the contact boundary. These surface 
discontinuities are associated with out-of-plane line forces. Such forces are not considered in 
the present framework. We thus consider a contact angle of 180°. Initially, prior to loading 
and contact, the droplet is spherical. We denote the initial radius R, and the initial volume 
Vo := 47ri?'^/3. The water inside the droplet is considered incompressible such that the vol- 
ume remains constant during deformation. The weight of the water causes hydrostatic pressure 



loading of the membrane leading to contact with the substrate. Fig. 11 shows the computed 



droplet deformation for various gravity values considering different EE formulations. As the 
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Figure 9: Cushion contact: configurations at pg = {0, 0.2, 0.4, 0.8}^/i?. The coloring shows 
the stress invariant /i = tr <t = cr? . 



gravity level increases the droplet spreads out on the substrate. Formulation (78) is used for 
the computations. Stability parameter p is set to 0.005 7. For this the errors in the surface 
tension are less than 2.5% in the case of quadratic Lagrange elements (Fig 11, top.) and less 
than 1% in the case of cubic T-splines (Fig [Tl| bottom). In both cases the errors increase 
along with pg. In the case of Lagrange elements the largest errors are found at the element 
boundaries, where the formulation is only C^-continous. For the T-spline case, which is 
continous over the entire surface (appart from two degenerate points at the top and bottom), 
the error is uniformly spread over the surface. It is remarked that the presented droplet model 
is much more general than the classical FE droplet formulation of Brown et al. (1980). This is 
discussed in detail in a forthcoming publication. 



5 Conclusion 



A novel computational formulation that is suitable for both solid and liquid, i.e. surface-tension- 
driven, membranes is presented. The theory, outlined in Sec. [2j is based on the differential 
geometry of curved surfaces, allowing for a very general formulation that accounts for large de- 
formations and general material laws. Curvilinear coordinates are used to formulate the surface 
geometry, kinematics, constitution, and balance laws. The governing strong and weak forms 
are split into the in-plane and out-of-plane parts, allowing the use of different approximation 
techniques for both parts and an elegant treatment of liquid membranes. Also, the consideration 
of deformation-dependent pressure loading comes naturally within the proposed formulation. 
Various constraints imposed upon the membranes can be handled by the theory, including vol- 
ume, area, and contact constraints. 

The membrane formulation is discretized using nonlinear finite elements. This results in a very 
efficient formulation that only uses three degrees of freedom per surface node and avoids the use 
of local cartesian coordinate systems and the transformation of derivatives. This is discussed 
in Sec. O 
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Figure 10: Growing liquid droplet: (a) pressure-volume relation for V £ [1/8 4]Vq (FE result 
for 12 quadratic FE with ^ = O.OI7); (b) convergence behavior for V = iVo and various fi, 
considering 3x3 Gaussian quadrature points. 



The capabilities of the formulation are demonstrated by several challenging examples in Sec. |4] 
Linear Lagrange, quadratic Lagrange, quadratic NURBS, and cubic T-spline finite elements 
are considered for the discretization. Constraints are imposed using the Lagrange multiplier 
method in the case of volume constraints and the penalty method for contact constraints. The 
inflation of a balloon and the growth of a droplet are used to validate the solid and liquid 
membrane formulations and they both yield excellent results. Comparing the different finite 
element types, the examples show that large accuracy gains lie between linear and quadratic 
Lagrange, and between quadratic Lagrange and isogeometric finite elements. 

The presented membrane formulation has been successfully applied to liquid droplets in this 



paper, but a rigorous analysis is still needed to assess the approach proposed in Eq. (78). 
Another important extension to the present formulation is the inclusion of bending stiffness, 
which can be present in both fluid and solid films. In the latter case this should lead naturally to 
a rotation-free shell formulation, which can be suitably handled by isogeometric finite elements. 
A further interesting extension is the consideration and development of different membrane 
material laws. Such a development is especially important for the case of biological membranes, 
which are often characterized by complex material behavior. 



A Consistent linearization of various quantities 

For Newton's method we need to linearize the kinematical quantities of the discrete system at 
X in the direction Ax. This is done at the FE level. 



A.l Linearization of 



According to Eq. (55) we have 



Aa„ = N,Q, Axe 



(80) 
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Figure 11: Liquid droplet contact: droplet deformation for pg = {1, 2, 4, 8, 20}^ /E?, each 
with identical volume and (i = 0.005 7: (a) quadratic Lagrange elements, (b) cubic T-spline 
elements. The color shows the membrane stress /i/2 normalized by the surface tension 7. In 
theory /i/(27) = 1. 



A. 2 Linearization of Oq,^ 

With definition ([s]) follows 



Acq^ = (tto, • N,/3 + ap ■ N,„) Axe 



(81) 



A. 3 Linearization of J 



The change A J can be written as 



where 



Thus 



A J = Attn 



daa 



dJ 
daa 



A J = Ja" • N,Q Axe 



(82) 

(83) 
(84) 



A. 4 Linearization of a"^ 

From Eq. ([5| and the formula 



a 



(85) 
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where 



is the unit alternator, we find 



1 
-1 



(86) 



(87) 



with 



A. 5 Linearization of nda 

The surface normal n appears together with the area element da and it is convenient to linearize 
them together. According to Eqs. ^ and (20) we have 



ndo = ai X a2 dD . 



Hence 



A(nda) = ^ i^N^i Axi x 02 + ai x N^2 Aa;/j dD . 

Expanding Axj into Axj = Axf + Ax" n we find 

Axj X a2 = Ja{nAx} — Axj) = Ja[n (>S) — (>Si n) Axj , 
ai X Ax I = Ja (n Ax'j — a? Axf) = Ja{n® a? — a? ® n) Axi , 



where J a = \/ det = da/dD. Thus 

A(n da) = [n ® a"' - a"" ® n) N_q, Axe da 

A. 6 Linearization of r"'^ 

For the solid model according to Eq. ( |32[ ) we have 

Ar'"/^ = ^T(2J-3 AJa"/^ - J-^ Aa"/^) , 
which can be rewritten into 



with 



Ar°^ = c^'^T^ • N 5 Axg , 



Note that the tensor [d^^'^l^], like [m"'^'^''], posses both major and minor symmetries. 
For the liquid model according to Eq. ([33]) we have 

Ar"^ = 7(AJa"'^ + J Aa"^) , 

which can also be written in the form ( |94[ ), where now 



(89) 
(90) 

(91) 

(92) 



(93) 

(94) 
(95) 



(96) 



7J -(e^^e'^'^ + e^'^e^^) 



(97) 
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B Finite element tangent matrices 



B.l Tangent matrix associated with 



int 



The internal force vector f , given in Eq. ( 63 ) , yields 



Aff^t = [ Ar"^ dAxe+ [ r"'^ N,;^ dA Ax^ . (98) 



In view of Eq. ( 94 ) , we can write 

Afi-„, = (kU + k^co)Axe (99) 
where we have introduced the material stiffness matrix 

k^at = / c"^^'^ N^, (a^ » a^) N,^ dA (lOO) 

and the geometric stiffness matrix 

k|eo= / ISi^T'^^f^^^dA (101) 

J fig 

Both these matrices are symmetric for the two constitutive models considered here. For those 
models, the terms in k^g^^ should be multiplied-out a-priory to obtain an efficient implemen- 
tation. If we consider splitting fj^^. into ^^^^ and fj^^.^ additional stiffness terms are picked up. 
These are reported in a forthcoming publication. 

B.2 Tangent matrix associated with Gl^^. 



From Eq. (67), for dead /q and t, we have 

Af|^t=/" N^nApda+ /" N^pA(nda) . (102) 



The first term is only required for hydrostatic loading according to Eq. (74). Here we find 



Ap = -pg'NAxe. (103) 
Contribution A(nda) is given by Eq. (92). As a result, 

k^^t = - / p'N^n(g>g'Nda+ p N'^ (n <^ a" - a" n) N,^. da . (104) 

B.3 Tangent contributions associated with the volume constraint 

If the volume constraint = O is active, we need to account for the unknown Lagrange multiplier 
Pv in the linearization. For the external forces we now have 



Af|,t = k|xtAxe + l|,tAp. , (105) 

with 
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llxt = ir^= / N^nda. (106) 



Further, at the element level, 
with 



= AXe 



(107) 



h^ = M = l/ n-Nda + l [ a;- (n®a"-a"®n)Nada . (108) 
The preceding contributions can be arranged into the elemental tangent matrix 



:= 



i^e _ T^e _ie 
*-int "-ext ^ext 



(109) 





which describes the change in fjnt — fext and gy due to changes in position Xg and pressure Pv 
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